R Notebook

Function to calculate the distribution function of the longest run and its mean. p - probability of success, n - number of trials.

library(plotly)
library(ggplot2)
longest_run_distribution <- function(p, n){
  mean = 0 
  E <- function(x){-(p^x) * exp(-(p^x)) * log(p, exp(1))}
  find_area <- function(l){
    lower_bound = l - log(n*(1-p), 1/p)
    upper_bound = lower_bound + 1
    area = integrate(E, lower_bound, upper_bound) 
    return(area[["value"]])
  }
  probabilities <- vector()
  labels <- vector()
  for(l in 1:n){
    area = find_area(l)
    
    # calculating mean
    mean = mean + area * l
    
    if(area > 0.001){
      probabilities <- c(probabilities, area)
      labels <- c(labels, toString(l))
    }
  }
  print("Mean:")
  print(mean)
  
  # plotting
  df <- data.frame(text=labels,
                  num=probabilities)
  
      plot_ly(data = df,
       x = ~text,
       y = ~num,
       type = "bar",
       color = I("darkred")
  ) %>%
  layout(
       title = "Distribution of the longest run lenghts",
       xaxis = list(title = "Length of the longest run",
       categoryorder = "array",
       categoryarray = ~text),
       yaxis = list(title = "Probability", range = c(0,max(probabilities)))
  )
}

Function to calculate the empirical distribution function of the longest run of Heads in n tosses of a fair coin and its mean.

longest_run_empirical_distribution <- function(n){
  trials = 10000
  longest_lens = vector()
  for (i in 1:trials) {
    coin <- sample(c("H", "T"), n, replace = TRUE)
    coin.rle <- rle(coin)
    longest_run_len = tapply(coin.rle$lengths, coin.rle$values, max)[1][["H"]]
    longest_lens <- c(longest_lens, longest_run_len)
  }
  
  #calculating mean
  mean = sum(longest_lens)/trials
  print("Mean:")
  print(mean)
  
  #plotting
  plot_ly(x=longest_lens, histfunc='sum', type = "histogram", color = I("orange")) %>%
  layout(title = "Empirical distribution of the longest run lenghts",
  xaxis = list(title = "Length of the longest run"),
  yaxis = list(title = "Frequency"))
  
}

Distribution of lenght of the longest run of Heads for n = 100 tosses of a fair coin:

p = 1/2
n = 100
longest_run_distribution(p,n)

[1] "Mean:"
[1] 5.976601

Empirical distribution of lenght of the longest run of Heads for n = 100 tosses of a fair coin:

longest_run_empirical_distribution(n)

[1] "Mean:"
[1] 5.9898

Distribution of lenght of the longest run of Heads for n = 1000 tosses of a fair coin:

p = 1/2
n = 1000
longest_run_distribution(p,n)

[1] "Mean:"
[1] 9.298531

Empirical distribution of lenght of the longest run of Heads for n = 1000 tosses of a fair coin:

longest_run_empirical_distribution(n)

[1] "Mean:"
[1] 9.288

Distribution of lenght of the longest run of Heads for n = 10000 tosses of a fair coin:

p = 1/2
n = 10000
longest_run_distribution(p,n)

[1] "Mean:"
[1] 12.62046

Empirical distribution of lenght of the longest run of Heads for n = 10000 tosses of a fair coin:

longest_run_empirical_distribution(n)

[1] "Mean:"
[1] 12.6057

Thus, we can see that cdf corresponds to the empirical cdf. An expected value obtained calculating numerically is close to the true mean. This is justified by the Law of large numbers. LLN says that the average of the results obtained from a large number of trials should be close to the expected value, and will tend to become closer as more trials are performed.